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Abstract. We present a computer simulation study of a (6,12)-Lennard-Jones fluid 
confined to a slit pore, formed by two uniform planes. These interact via a (3,9)- 
Lennard-Jones potential with the fluid particles. When the fluid approaches the 
liquid-to-solid transition we first observe layering parallel to the walls. In order to 
investigate the nature of the freezing transition we performed a detailed analysis of the 
bond-orientational order parameter in the layers. We found no signs of hexatic order 
which would indicate a melting scenario of the Kosterlitz-Thouless type. An analysis 
of the mean-square displacement shows that the particles can easily move between the 
layers, making the crystallization a 3d-like process. This is consistent with the fact 
that we observe a considerable hysteresis in the heating-freezing curves, showing that 
the crystallization transition proceeds as an activated process. 

PACS numbers: 



Submitted to: J. Phys.: Condens. Matter 



2D versus 3D Freezing of a LJ Fluid in a Slit Pore 



2 



1. Introduction 

Understanding tlie structure and dynamics of confined fluids is important for processes 
such as wetting, coating, and nucleation. The properties of a fluid confined in a pore 
differ significantly from the bulk fluid due to finite size effects, surface forces and reduced 
dimensionality. In this work we report a study of one of the simplest models that is still 
capable of reproducing the thermodynamic behavior of classical fluids, the Lennard- 
Jones (LJ) system. The LJ potential is an important model for exploring the behavior 
of simple fluids and has been used to study homogeneous vapor-liquid, liquid-liquid and 
liquid-solid equilibria, melting and freezing. It has also been used as a reference fluid 
for complex systems like colloidal and polymeric systems. 

The vapor-to-liquid transition in confined systems has been studied intensively, 
and it is well understood (see [1] and references therein). In this article we will discuss 
the liquid-to-solid transition in a slit pore and the process of the development of the 
solid phase. In the liquid phase, confinement to a slit induces layering at the walls. 
One could imagine this effect to facilitate crystallization. And indeed it is known that 
depending on the strength of the particle-wall interaction different scenarios of freezing 
exist [21 [3]. If the walls are strongly attractive, crystallization starts from the walls 
and at a temperature higher than without confinement. If, however, the walls are 
strongly repulsive, crystallization starts from the bulk at a temperature lower than 
without confinement. A well-distinguished layer of particles at the wall can also, to 
some extent, be treated as a 2d system. This raises the question, whether freezing 
of such a layer proceeds via the Kosterlitz-Thouless-Halperin-Nelson- Young (KTHNY) 
mechanism |H O El [7], meaning that the liquid turns into a crystal going through a 
hexatic phase This question has been studied for rather narrow pores (up to 7.5 
diameters of a fluid particle). It was found that a hexatic phase exists between liquid 
and crystal only in the contact layers at the walls [9]. As a consequence, in a pore that 
can accommodate only a single layer the transition is of the KTHNY type. However, 
with increasing width the behavior changes to a first order transition jlOj . 

Here we investigate an attractive pore that is significantly wider, namely 20 
diameters of a fluid particle. Studying the bond-orientational order parameter within 
the layers we observe no sign of a hexatic phase. An analysis of the mean- 
square displacement shows that the particle diffuse between the layers. Hence, the 
crystallization proceeds as a 3d process, as is also suggested by the noticable hysteresis 
loop in the heating-freezing curve. 

The article is structured as follows. First we describe our simulation method. In 
section [3l we present the results, followed by a discussion, and in section HI we conclude 
with a summary of the presented results. 
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2. Simulation method 

We performed molecular dynamics (MD) simulations in the isothermal ensemble (NVT), 
i.e. the number of particles N, the volume V and the temperature T were fixed. The 
system consists of 8000 particles confined between two structureless walls. The particles 
interact via the LJ-potential 



The interaction between walls and particles is characterized by a LJ-potential integrated 
over semi-space: 



The particle-particle interaction was cut off at a distance = 2.5a and the wall-particle 
interaction at a distance = 4. Oct (the wall-particle potential is wider and deeper then 
the particle-particle potential). Using for the wall-particles interaction a Steelle potential 
[llj or just a (4,10)-LJ potential does not influence the results qualitatively. For the 



following we will use e as the unit of energy, a as the unit of length and r = y 1 ■ cr^/e as 
unit of time (i.e. use the particle mass as the unit of mass); consequently, temperatures 
are given in multiples of e/k-QT. The simulations were performed in a cubic box with 
periodic boundary conditions in the x and y direction. The walls were positioned at 
z = and z = = 20a. The size of the simulation box was = Ly = = 20a. 
We used standard Nose-Hoover and Langevin thermostats to keep the temperature 
constant [121 [131 ^M- For the Langevin thermostat, the friction was always chosen as 
r = mT~^ = 1 [13]. For the Nose- Hoover thermostat, we set the effective mass Mg = 0.5. 
We simulated a cooling curve starting out from a random configuration at T = 3.0 and a 
melting curve starting out from a face-centered-cubic configuration at T = 1.2. Far away 
from the transition, the temperature was changed by AT = 0.1 from one simulation run 
to the next, while close to transition we used a smaller increment /decrement, AT = 0.01. 

The simulations were performed with a timestep of At = 0.005 and let run for 
1.0 X 10^ MD steps for equilibration and for 2.5 x 10^ for sampling. In order to compute 
the mean square displacement of the particles, a smaller timestep At = 0.002 was chosen, 
and to avoid influence of thermostat on the dynamics of the system we switched to the 
NVE ensemble after equilibration (keeping the total energy of the system E constant). 
We monitored the temperature, which fluctuated around a mean value practically equal 
to the temperature T in the NVT ensemble. Pressures were obtained from the virial 
expansion [15] omitting corrections for the cut-off in the potential. For parts of our 
simulations the software package ESPResSo version 2.04s was used [T6] . 

3. Results and Discussion 




(1) 




(2) 




In Figure [H the pressure-temperature curves for heating and cooling are shown. There is 
a considerable hysteresis, which indicates that the system has to overcome a free energy 
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barrier when transforming from one phase to the other. Several runs were performed 
both for heating and for coohng. The coohng hnes coincide, whereas the temperature at 
which the melting process starts fluctuates. In the Figure [1] we show the outer borders 
of the hysteresis region. 

As our system has attractive walls, crystallization should start from the walls [2]. 
This can be clearly seen in the snapshot (Figure [2]) that was taken 90 MD steps after 
equilibration had started. Only another 90 MD steps later the system completely 
crystallized. One can also see that no crystallization process has started at the right 
wall yet, demonstrating that this event is an activated process. 

As it was shown in [2] the width of hysteresis depends on the distance between the 
layers. If the distance differs considerably from the lattice constant of an ideal L J crystal 
(0.916), then the hysteresis will be more pronounced. For our system the distance is 
0.87 in the bulk, and correspondingly the hysteresis is quite wide. 

In order to investigate the phase transformation process, we now turn to the effects 
the walls have on the structure of the fluid. Figure [3] shows number density profiles g{z) 
for Lz = 20. Oct in the liquid and the solid phase. In the liquid phase, the maxima of the 
peaks follow an exponential law A [exp{—Bx) + exp(— (L^ — x)B)] + Pmid, where pmid is 
the density in the middle of the box. Figure H] shows the behavior of the coefficient B 
with temperature. It can be seen that the values of B decrease more or less linearly at 
first, i. e. the number of layers increases and they become more pronounced. As soon 
as we enter the regime of the hysteresis at T = 2.0, B becomes almost constant (within 
the error of the simulations). This shows that the structure of the density profile does 
not change, no new layers appear and the system is trapped in the undercooled state. 

As the liquid forms layers, one could assume that the transformation proceeds inside 
the layers via a KTHNY transition. In order to test this assumption, we now turn to 
the structure within the layers: To characterize the transitional order in one layer in 
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Figure 2. Snapshot of the system in the early stage of crystaUization at T=1.60. In 
this specific example it crystallizes from the left wall. Walls are not represented, but 
are to the right and to the left of the box. The part on the left is already a crystal 
while the right side is still disordered. 
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Figure 3. Density profile g{z) of the liquid phase at T=1.70. The heights of the 
peaks are fitted by an exponential function. The inset shows g{z) for the solid phase 
at T=1.60. 



2D, we calculate the pair correlation function: 

9{r) = g-'{J:Sir^nr,-r)) (3) 

where g is the number density of particles in each layer. In Figure [5] the 2D radial 
distribution functions for the first and third layer (seen from the wall), the bulk part of 
the liquid and the first three layers of the solid phase are shown. The structure within 
the layers of the liquid becomes less pronounced as we move further away from the walls 
and is barely visible in the center of the box. 
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Figure 5. 2D pair correlation function in the first and third layer and the bulk part 
of the system in the the liquid state at T = 1.65. The inset shows g{r) for the first 
three layers at the wall for the solid state at T = 1.60 



Next we consider the bond-orientational order [8j: we define the local bond- 
orientational order parameter of particle j in layer m at a position Xj as 

C(x,) = ^Ee^''^'= (4) 

i k=l 

where Nj is the number of neighbors of particle j within layer m, the sum is over the 
neighbors k of j within m, and 6jj^ is the angle between an arbitrary fixed axis and 
the line connecting particles j and k. The order of the m-th layer \Ef^ is defined as the 
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Figure 6. 2D bond-orientation order parameter for liquid states depending on 
the layer m for temperatures 1.70, 1.80, 2.00, 2.20, 2.40 (from top to bottom). The 
inset shows ^™ for solid states for temperatures 0.80, 1.00, 1.30, 1.50, 1.60 (from top 
to bottom). 
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Figure 7. 2D bond-orientation order parameter of the first layer from the wall as 
a function of the temperature. 



average over ?/^^(xj) for all Nra particles within the layer 

1 Nrr, 

^r = ^IEC(x,)l . (5) 

Figure [6] shows \E'^ for various temperatures. When approaching the transition, the 
bond-orientational order close to the wall increases. 

The temperature dependence of \l/g for the first layer of particles at the wall is 
shown in Figure [71 It clearly "jumps" i. e. is discontinuous at the transition. 

If the crystallization proceeded purely within the two-dimensional layers, one would 
observe a hexatic phase, which is characterized by a power-law decay of the correlation 
of the bond-orientational order 



(6) 
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Figure 8. Bond-orientational correlation function for the liquid state just before 
freezing at T = 1.67 for the layer closest to the wall. The solid line is the result of 
an exponential fit. There is no signature of a hexatic phase. The inset shows gg for 
the crystallin state at the next available lower temperature T = 1.66, right after the 
transition. 

where the average is taken over all particles within a layer whose positions x are a 
distance r apart. 

Figure [8] shows g^ir) for the first layer at T = 1.67 and T = 1.66. The system 
jumps from the 2d-liquid phase into the 2d-solid without visiting a hexatic phase first. 
Hence the crystallization process is not of the KTHNY-kind. 

To find out why the crystallization process is 3d-like despite the layering, we now 
consider the particles' dynamics. One of the obvious characterictics is to estimate how 
long particles on average stay in the layer closest to the wall. The easiest way to 
estimate this is to calculate how many particles of those which were in the layer at time 
remained there at the time t. From Figure [9] we can see that the ratio of particles that 
remain in the layer decreases exponentially with time. Fitting it with exp(— t/r) we 
obtain the average lifetime r of a particle in the layer (Figure fTOi) . It increases linearly 
with the decrease of temperature and is then fluctuating around the mean value in the 
hysteresis region. As we observed already for the density, the behavior of the system in 
the hysteresis region does not change much during cooling. 

To characterize the mobility of the particles we calculated the mean square 
displacement (MSD). As the system forms layers, we calculate the MSD parallel and 
perpendicular to the wall separately. Looking at the plane parallel to the wall (Figure [TTl) 
while approaching crystallization, we observe that the particles in the layer closest 
to the wall are a little faster than the particles in the bulk, despite the fact that 
the crystallization typically starts from here. We take this as another hint that the 
crystallization proceeds as a 3d-process, and does not first start within the layer closest 
to the wall. The behavior of the particles does not change significantly on approach of 
the crystallization as they enter the metastable region. 

The mean square displacement measured perpendicular to the wall (Figure [11]) 
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Figure 9. Ratio of particles in the layer closest to the wall that stayed there from 
time (iV(0)) until time t {N{t)) at different temperatures above the phase transition. 




Figure 10. Average lifetime of particles in the outmost layer as a function of 
temperature. In the hysteresis region the lifetime does not change significantly. 



shows that after the balhstic regime for a while particles are trapped in the layer and 
then start leaving it. It is not meaningful to calculate the diffusion coefficient in our 
system, because the particles do not stay long enough in a layer for the MSD to enter 
the linear regime. 



4. Conclusions 



We reported on a molecular dynamics study of the liquid-to-solid transformation of a L J 
fluid in a wide slit pore. Although the confinement induces layering in the liquid phase 
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Figure 11. Mean square displacement for T = 1.64. The MSD parallel to the wall is 
almost identical in the layer and in the bulk with the particles in the layer even being 
slightly faster. The MSD perpendicular to the wall show a clear trapping effect. 



close to the walls, we do not find a successive, layerwise crystallization. Crystallization is 
still a 3d process, and, in particular, no hexatic phase was observed in the layers closest 
to the wall, excluding the possibility of a 2D KTHNY-like crystallization within the 
layers; in fact, the mobility of particles in the layers is higher than their mobility in the 
bulk. Nevertheless, we find that crystallization in the system practically always starts 
from the walls, i. e., the walls facilitate crystallization. And although crystallization 
is an activated process similar to 3d crystallization, we observe a smaller hysteresis, 
indicating a reduced nucleation barrier as compared to bulk crystallization. 

AUtogether, our simulations suggest that the nucleation of the LJ fluid close to a 
planar wall does not significantly differ from the nucleation in the bulk, although with 
a smaller nucleation barrier. This can however be easily understood as an effect of the 
strongly increased density in the layers close to the confinement. 
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